%计算单个电机在某一损耗一段时间下的温升情况；
%输入input1: [定、转子损耗]kW, input2: [定、转子温度]℃, input3: 海拔m, input4: 时间s
%在原函数基础上将环境温度修改为输入参数
function [theta_s,theta_r] = lptn2nodes_sr5(PLoss, Temp, h, time, envTheta)
% 热路参数由GA辨识所得
Rs = 25.21994135;      % 定子热阻
Cs = 122.77;      % 定子热容
R_sr = 23.998;    % 气隙热阻
Rr = 87.781;      % 转子热阻
Cr = 89.345;       % 转子热容

st_ploss = PLoss(1);
rt_ploss = PLoss(2);

st_temp = Temp(1);
rt_temp = Temp(2);

% envTheta = 20-6*h/1000;    % 根据海拔计算环境温度
matC = [-1/Cs*(1/Rs+1/R_sr) 1/(Cs*R_sr);...
          1/(Cr*R_sr)   -1/Cr*(1/Rr+1/R_sr)];     % 移项后矩阵
matP = [1/Cs  0;...
         0  1/Cr];   
matE = [1/(Cs*Rs); 1/(Cr*Rr)] * envTheta;
K = 1;  % 散热风量系数 standard=1

%计算冷却功率 二阶多项式近似Motor-CAD结果 负值 散热考虑为负功率
rotorSubEnv = rt_temp - envTheta;    % 计算风扇对转子冷却功率
if rotorSubEnv > 5
    rotorWindLoss = K * (1439 +  -0.4851*h + (-61.3)*rotorSubEnv + 0.02029*h*rotorSubEnv + ...
        -0.4169*rotorSubEnv^2 + (-4.255e-05)*h*rotorSubEnv^2 + 0.001287*rotorSubEnv^3) / 1000;
else
    rotorWindLoss = 0;
end
statorSubEnv = st_temp - envTheta;   % 计算风扇对定子冷却功率
if statorSubEnv > 5
    satorWindLoss = K * (913.7 + (-0.3241)*h + (-72.66)*statorSubEnv + ...
        0.01186*h*statorSubEnv + (-0.1177)*statorSubEnv^2) / 1000;
else
    satorWindLoss = 0;
end
%热耗矩阵
loss_mat = [st_ploss + satorWindLoss;...
            rt_ploss + rotorWindLoss];
% 瞬态热流微分方程矩阵形式
[T,Y] = ode45(@(t,y) matC*y + matE + matP*loss_mat, [0,time], Temp);

endTemp = Y(end,:);
theta_s = endTemp(1);
theta_r = endTemp(2);